F(* Mathematica Package *)

BeginPackage["InterfacialThinFilmNLS`", {"TransferMatrixFormalism`","KramersKronigOscillators`"}];

SetDirectory[DirectoryName @ $InputFileName];

Unprotect["`*"]; 
ClearAll["`*"];

PhaseMatchAngle::usage = "imported"
Incidence::usage = "imported";
InterfaceSide::usage = "imported";
AzimuthAngle::usage = "imported";
TransferCoefficientMethod::usage = "imported";
LocalFieldsFromLimits::usage = "imported";
SheetThickness::usage = "imported";
DetectionSide::usage = "imported";
FresnelL::usage = "imported";
NonlinearPolarization::usage = "imported";
InputTransferC::usage = "imported";
OutputTransferC::usage = "imported";
GenerateInterfacialModel::usage = "";

Get[FileNameJoin[{"Usages", "InterfacialThinFilmNLS.m"}]]

Begin["`Private`"]

<< ThinFilmNLSupdater`
CheckUpdate[DirectoryName @ $InputFileName]

Unprotect["`*"]; 
ClearAll["`*"];

PhaseMatchAngle[OmegaTheta:{{_,_}..}] := ArcSin[Total[OmegaTheta[[All, 1]] * Sin[OmegaTheta[[All, 2]]]] / (Total @ OmegaTheta[[All, 1]])]
PhaseMatchAngle[OmegaThetaIndex:{{_,_,_}..}, IndexWM_]:=ArcSin[Total[Times /@ OmegaThetaIndex[[All, {1, 3}]] * OmegaThetaIndex[[All, 3]] * Sin[OmegaThetaIndex[[All, 2]]]] / (IndexWM * Total@OmegaThetaIndex[[All,1]])]
(*PhaseMatchAngle[OmegaThetaIndexPhi:{{_,_,_,_}..}]= *)
PhaseMatchAngle[___]/;Message[General::badargs,PhaseMatchAngle]:="unevaluated"

NonlinearPolarization::baddims = "Product of last two dimensions of \[Chi] = `1` should match product of lengths of all input fields = `2`."
NonlinearPolarization[chi_,Ein:{{_,_,_}..}]/; Times @@ Rest[Dimensions[chi]] == Times @@ (Length /@ Ein) := 
	Total[#, Infinity] & /@ (Times[#, Fold[Outer[Times, #1, #2]&, First @ Ein, Rest @ Ein]]& /@ chi)
NonlinearPolarization[chi_,Ein:{{_,_,_}..}]/;Message[NonlinearPolarization::baddims, Times @@ Rest[Dimensions[chi]],Times @@ (Length /@ Ein)]:= "not evaluated"

	
(*Input Field Transfer Coefficient*)
Pr = {{Cos[#], 0, 0},
	  {0, 1, 0},
	  {Sin[#], 0, 0}}&;(*[Theta_]*)

J = DiagonalMatrix[{1, 1, (#1/#2)^2}]&;(*[n_,nbar_]*)

Zrot[AzimuthAngel_]:= RotationMatrix[AzimuthAngel, {0, 0, 1}];(*[phi_]*)
Zrot[0] = IdentityMatrix[3];


Get[FileNameJoin[{"DownValues", "InputTransferC.m"}]];

(*Output Field Transfer Coefficient*)
(*need rs and rp*)
rs = Evaluate @ FresnelA["rs",#1, #2, #3, #4]&;
rp = Evaluate @ FresnelA["rp",#1, #2, #3, #4]&;

pj = (I #1 Sec[#2]) / #3 &;(*[omega, theta, n]*)

f = DiagonalMatrix[{rp[#1, #2, #3, #4], rs[#1, #2, #3, #4], 0}]&;(*[thetai, thetaj, ni, nj]*)

FresnelL["v+", {thetavm1_, thetav_}, {nvm1_, nv_}, nbarv_] := 
	(Transpose[Pr[thetavm1]] + f[thetavm1, thetav, nvm1, nv] . Transpose[Pr[-thetavm1]]) . J[nv, nbarv];

FresnelL["v-'", {thetavm1_, thetav_}, {nvm1_, nv_}, nbarv_] := 
	(Transpose[Pr[-thetav]] + f[thetav, thetavm1, nv, nvm1].Transpose[Pr[thetav]]).J[nvm1, nbarv];
	
FresnelL["v+", thetavm1_, thetav_, nvm1_, nv_, nbarv_] := FresnelL["v+", {thetavm1, thetav}, {nvm1, nv}, nbarv];
FresnelL["v-'", thetavm1_, thetav_, nvm1_, nv_, nbarv_] := FresnelL["v-'", {thetavm1, thetav}, {nvm1, nv}, nbarv];
FresnelL["v+", thetavm1_, nvm1_, nv_, nbarv_] := FresnelL["v+", {thetavm1, RefractionAngle[thetavm1,nvm1, nv]}, {nvm1, nv}, nbarv];
FresnelL["v-'", thetavm1_, nvm1_, nv_, nbarv_] := FresnelL["v-'", {thetavm1, RefractionAngle[thetavm1,nvm1, nv]}, {nvm1, nv}, nbarv];
FresnelL[___]/;Message[General::badargs,FresnelL]:= "unevaluated";



Get[FileNameJoin[{"DownValues", "OutputTransferC.m"}]];

Get[FileNameJoin[{"DownValues", "GenerateInterfacialModel.m"}]];


Protect["`*"]; 

End[]

Protect["`*"];

ResetDirectory[];

EndPackage[];